load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  expt  = (/"hist-aer","hist-GHG","hist-nat","historical"/)

  model = (/"ACCESS-CM2","ACCESS-ESM1-5","BCC-CSM2-MR","CNRM-CM6-1","GFDL-ESM4",\
            "HadGEM3-GC31-LL","IPSL-CM6A-LR","MIROC6","MRI-ESM2-0","NorESM2-LM"/);10 CMIP6 models
 
  fread = addfile("data/pr_eof_1951-2020_JJAS_GPCC.nc","r")
  eof0  = fread->eof_regress(0,:,:)
  eof0  = eof0*(-1.)

  fread  = addfile("/WORK/jiangj/data/obs/rainfall/GPCC/full_data_monthly_v2022_10.nc","r")
  lat    = fread->lat({25:42})
  lon    = fread->lon({65:105})

  eof = linint2_Wrap(eof0&lon,eof0&lat,eof0,False,lon,lat,0); regrid to 1°×1°

  fread1 = addfile("/WORK/jiangj/data/GMTED2010/GMTED2010_15n030_0125deg.nc","r")
  h      = short2flt(fread1->elevation)*1.
  lat0   = fread1->latitude
  lon0   = fread1->longitude

  h!0    = "lat"
  h!1    = "lon"
  h&lat  = lat0
  h&lon  = lon0
  h      = h/1000.

  hh  = h({25:42},{65:105})
  h  := linint2_Wrap(hh&lon,hh&lat,hh,False,lon,lat,0)
  eof = where(h.ge.2.5,eof,eof@_FillValue); only keep the pattern over HMA

  do i=0,9
    print(model(i))

    do j=0,3
        print(expt(j))
        
        dilo = "data/CMIP6_DA/"+model(i)+"/"+expt(j)+"/mon/pr/"
        fils = systemfunc("ls "+dilo+"pr_Amon_"+model(i)+"_"+expt(j)+"*185001-202012.nc")
        nf   = dimsizes(fils)

        ts = new((/nf,70/),"float")

        do m=0,nf-1
          fread = addfile(fils(m),"r")
          time  = cd_calendar(fread->time,1)
          tStrt = ind(time.eq.195101)
          tLast = ind(time.eq.202012)      
          pr0   = fread->pr(tStrt:tLast,{20:50},{60:110})
          pr0   = pr0*86400*30.; mm/month
          pr1   = linint2_Wrap(pr0&lon,pr0&lat,pr0,False,lon,lat,0)
          pr    = month_to_season(pr1,"JJA")
          pr    = (pr1(5::12,:,:)+pr1(6::12,:,:)+pr1(7::12,:,:)+pr1(8::12,:,:))/4.; summer mean
          pr    = where(ismissing(conform(pr,eof,(/1,2/))),pr@_FillValue,pr); only keep the precipitation over HMA

          ts(m,:) = regCoef_n(eof,pr,(/0,1/),(/1,2/));timeseries
        end do 

        tm = dim_avg_n_Wrap(ts,0) ; the ensemble mean of all members 

        foutname = "data/pr_pc1_1951-2020_JJAS_"+expt(j)+"_"+model(i)+".nc"; timeseries for each model
        system("rm -rf "+foutname)
        fout = addfile(foutname,"c")
        fout->ts = ts
        fout->tm = tm

        delete([/fread,pr0,pr1,pr,ts/])

        delete(fils)
    end do 

  end do 

end 